function [L,U,P] = gaussp(A)

% ELIMINAZIONE DI GAUSS CON PIVOT PARZIALE
% Dati di input:
%   A = matrice quadrata
% Dati di output:
%   L = matrice triangolare inferiore
%   U = matrice triangolare superiore
%   P = matrice di permutazione
% tali che PA = LU

[n,m] = size(A);
if n ~= m
   disp('ERRORE: la matrice di input non e'' quadrata')
   break
end

P = eye(n);

for k = 1:n-1

% ricerca del pivot e permutazione:
% si scambiano le righe k e rp di A e di P
   [piv,indm] = max(abs(A(k:n,k)));
   indm = indm + k - 1;
   A([k indm],:)=A([indm,k],:)
   P([k indm],:)=P([indm,k],:)
  

% controllo sulla grandezza dell'elemento pivotale:
if abs(A(k,k)) <= realmin
      disp('Un elemento pivotale e'' nullo')
      error('La matrice in input e'' singolare')
   end

   A(k+1:n,k) = A(k+1:n,k)/A(k,k);
   A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
end

L = tril(A,-1) + eye(n);
U = triu(A);
